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ABSTRACT 

There is growing evidence that the active galactic nuclei (AGN) associated with the 
central elliptical galaxy in clusters of galaxies are playing an important role in the 
evolution of the intracluster medium (ICM) and clusters themselves. We use high 
resolution three-dimensional simulations to study the interaction of the cavities cre- 
ated by AGN outflows (bubbles) with the ambient ICM. The gravitational poten- 
tial of the cluster is modelled using the observed temperature and density profiles 
of the Virgo cluster. We demonstrate the importance of the hydrodynamical Kutta- 
Zhukovsky forces associated with the vortex ring structure of the bubbles, and discuss 
possible effects of diffusive processes on their evolution. 

Key v^rords: galaxies: cooling flows - galaxies: nuclei - galaxies: active - galaxies: 
clusters: general - galaxies: clusters: individual: Virgo - methods: numerical 



1 INTRODUCTION 

It has been more then 40 years since the first detection of x- 
ray emission from a cluster of galaxies was ma de. The emis- 
sion w as detected from around the M87 galaxy (iBvram et al.l 
1 19661 ) at the centre of the Virgo cluster, whose proximity 
helped to get high resolution maps of the x - ray emitting gas 
llMatsushita et all [200^ : lYoung et aLll2002l : iGhizzardi et all 
I2OO4I . hereafter G04). Later observations revealed that many 
clusters are bright x-ray sources, with luminosities in the 
range lO*^" *^ ergs~^. The properties of the x-ray emission 
from clusters of galaxies were found to be most consisten t 
with thermal bremsstrahlung from hot gas l|Sarazinlll98^ . 
This implies that the observed ICM has an electron density, 
rio, in the range of lO"* ""'^ cm~^, and a temperature, Ta, 
of the order lO'^- K. 

The radiative cooling time of the gas in a cluster core 
due to bremsstrahlung can be as short as 10^ yr. Unless the 
gas is thermally supported, the cooling leads inexorably to 
an inflow of the cold gas onto the central galaxy (for a re- 
view of cooling flow theory see, e.(;.,lFabian 1994). However, 
many observations have since demonstrated both the lack 
of the cold gas deposits, and spectroscopically determined 
mass deposition rates up to an order of magnitude smaller 
than the rates predicted by the classical cooling flow model 
ijVoigt fc Fabia n 2004). 

This apparent contradiction has led to the sugges- 
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tion that the gas in the cluster cores is reheated. A 
number of different mechanisms were proposed to re- 
duce the cooling flow (or prevent it from forming in 
the flrst place): currently the most popular candidates 
are outflows from active galactic nuclei (AGN) (see, e.g., 
Churazov eral]|200ll : iBriiggen fc Kaised[200^ : iBirzan et all 
2004h : thermal conduction of heat from t h e outer re- 
' MedvedevI I2OOII: IVoigt et al.1 

looi); 



gions (see, e.g.. 



Narayan & 
2004; 



l2002l : IVoigt fc Fabian 2004; Pope ct al.1 120061 ): disruption 
and heatin g of the flow by sound waves and turbulence 
(see, e.g., iRuszkowski et al.1 |2004| ; iFuiita fc Suzukil l2005l : 
iFuiita et al.ll2004l): heat i ng by supernovae explosions (see. 



Voit fc Brvanll200ll : iPipino et al1l2002l ; iTang fc Wan3 

l2005h . Clear discrimination between the different models 
based on the observational data is difScult. Morphological 
features of the x-ray emission were used as the evidence in 



favour of a particular mechanism I Ruszkowski et al]|2004l ; 
iFabian et alll2003l : iBriiggen et"al] |2005l). The morphological 
resemblance is important, and although it is impossible to 
base a prove in favour of a particular mechanism on the 
morphology alone, it is important to understand how differ- 
ent physical processes affect the morphology of the observed 
features. 

Despite some recent progress in the underst a nding 
of the physical state of the ICM l|Enfihn et al.1 I2OO5I : 
lEnfilin fc Vo"gtll2006l : iLazarianllioO^ ). many key parameters 
remain not well known. Since laboratory experiments on 
plasma with conditions cl ose to those found in astronom- 
ical objects are still rare IjKeenan fc Rosdl2004h , all esti- 
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mates of such physical properties of the ICM as the val- 
ues of thermal diffusivity and viscosity are generally based 
on mod ification s of th e values calculated in the classical 
work of ISpitzeJ (|l962h . It is accepted that magnetic fields 
decrease the diffusivity, since they suppress movements of 
charged particles across the field lines. The same is be- 
lieved to be true for the value of the viscosity as well. Usu- 
ally a fraction of the Spitzer va lue is used to approximate 
thermal difi'usivity of the IC M (|Chandran fc Cowley 1998.: 
iNaravan fc Medvedevll200ll '). which assumes chaotic orien- 
tation of the magnetic field lines [i.e., magnetic fields lines 
are bent on scales smaller than the mean free path of an 
electron in the plasma). The value of the thermal conduc- 
tivity in turbulent ICM can be an or der of magn i tude l arger 
than the Spitze r value, as shown by ICho et al] (|2003l ) and 
iLazarianl (I2OO6I '). The exact extent of this amplification de- 
pends on the state of the turbulence in the ICM, which is 
difficult to estimate from the present observational data. 

The kinetic effects are likely to be important in the 
context of the ICM physics. CoUisionless interactions on the 
scale smaller than the mean free path can affect th e overall 
dynamics o f the AGN-ICM interaction (|Schekochihin et al.l 
I2OO6I . I2OO5I ). 

1.1 The Model 

While a model that accounts for the above mentioned phys- 
ical properties of the ICM is clearly necessary for a detailed 
understanding of its dynamics, we will use a simplified ap- 
proach in order to investigate certain aspects of the dynam- 
ics related to AGN activity. In this work we use a purely 
hydrodynamical model of the ICM, which includes radiative 
cooling to simulate conditions preexisting to AGN activity 
due to the cooling flow. 

We use numerical simulations to study the morphology 
of the AGN-blown bubbles. We do not simulate inflation of 
the bubbles by AGN outflows, instead we introduce them as 
perturbations of the temperature and density fields of the 
cluster at a certain point during the simulation (see >i2.3l 
for more details) . While the inflation of the bubbles can in- 
fluence their subsequent evolution, we are more interested 
in the time scales that are larger than the inflation time, 
and the effects of numerical techniques on the hydrodynam- 
ics of the bubbles. Although the effect of the velocity held 
produced by the jet inflation of the bubbles is likely to be 
important especially during the early stages of the evolution, 
here we study the influence of the velocity field induced by 
the buoyantly rising bubble and the fluid instabilities devel- 
oping on its surface on the morphology of the bubbles. The 
combined effect of the jet injection and buoyancy will be 
addressed in the future work. 

The gravitational potential of the cluster is modelled us- 
ing observed temperature and density profiles for the Virgo 
cluster as determined by G04, using XMM, BeppoSAX and 
Chandra data (see i ]2.2l for more details). The observational 
data allows for more accurate modelling of the gravita- 
tional potential the n a fitted theoretical profile {e.g., NFW, 
iNavarro et al.lll997h , as the gravity of the central galaxy can 
become dominant at the cluster core. 

The outline of the current work is as follows: in § [2] we 
discuss numerical framework and initial conditions for our 
simulations, in § [3] we present analysis of the simulations. 



In §|4]we discuss implications and limitations of the current 
work. Finally, we present our main conclusions in § [5] 



2 NUMERICAL FRAMEWORK 

For our numerical experiment we use FLASH (version 2.3) - 
a modular, adaptive-mesh (AMR), parallel simulation code, 
capable of handling general compressible fl ow problems 
found in many astrophysical environments IjFrvxell et al.l 
l200d ). 

The HYDRO module of the FLASH code solves Euler's 
equations in three dimensions. In the conservative form the 



equations are given by, 

g + V ■ (pv) = 0, (1) 

^ + V ■ (pv ® v) + VP = pg (2) 
at 

M + v((p£ + P)v)=p(v.g) (3) 



where p is the fluid density, v is the fluid velocity, P is 
the pressure, £ is the sum of the internal energy E, and the 
kinetic energy per unit mass, £ — i5+|vp/2, A is the cooling 
function (see §[2TT}, and (v ® v)ij = ViVj denotes a tensor 
product. 

The physical size of the computational grid is 10^* cm ~ 
322 kpc in each direction. Numerically the grid is con- 
structed using nested blocks of 16^ cells, with the ratio 
of sizes between the neighbouring blocks being either 1:1 
or 1:2. The minimum refinement level (minimum level of 
nested blocks with 1:2 ratio) is set to 3, which results in 
the largest cell size of 10^* cm/(2^~^16) = 10^* cm/64 « 
1.56 X 10^^ cm « 5 kpc. The maximum refinement level was 
set to 8, giving a minimum cell size of 10^* cm/(2*~^ 16) = 
10^^ cm/2048 f» 4.88 x 10^° cm « 0.16 kpc. 

2.1 Radiative Cooling 

We have developed a new module for FLASH, which imple- 
ments the cooling function, A (see Fig. [ij , to account for the 
radiative losses from the fully ionised plasma in the wide 
range of tem peratures, 4 < logT < 8.5. using the values 
tabulated bv [Sutherland fc Dopital (|l993l ) (metallicity was 
taken to equal half solar [Fe/H]= —0.5). For temperatures 
exceeding logT = 8.5 we calculate A using the formula for 
thermal bremsstrahlung (|McKee fc Cowi3ll977l ). 

E = --A, (4) 
P 

A = ncnpAw, (5) 

An = 2.5 X 10"^^r°-^ [ergs"^cm% (6) 

where E is the rate of the energy loss per unit mass, rio, rip 
are the electron and proton number densities, which for fully 
ionised gas with primordial abundances (mass of Helium is 
0.25 of total mass of the gas) relates to the mass density of 
the ICM as, 

ne = 1.167np, (7) 
p = 1.143nca, (8) 

where a = 1.661 x 10"^'* g is the atomic mass unit (~ rUp). 
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Figure 1. Dependence of the cooling function, An, on tempera- 
ture (see ? I2.1|I : thick Une — the current model; dashed line — the 
original FLASH cooling function; dotted line — the cooling function 
from McKee & Cowic (1977). 



The new cooling function accounts for line cool- 
ing in much greater detail, but it does not dif- 
fer significantly from the default FLASH module 
(source_terms/cool/radloss), see Fig. [1] The only 
serious departure from the default approximation occurs 
around temperature of T ~ lO'^ K, were the original module 
overestimates the energy loss by a factor of a few (~ 5). It 
is important to note, that in our simulations most of the 
gas has temperatures within this range, and therefore this 
difference is non-negligible. The radiative losses were not 
applied to the gas with temperatures below T = 10* K 
because the present scheme does not thermally resolve gas 
at such temperatures (cooling length becomes smaller than 
the numerical resolution of our grid). 

2.2 Initial Conditions 

We model the cluster gravitational potential using the data 
from G04. Using the de-projected temperature and electron 
number density profiles of the Virgo cluster. 



T{x) = Ti - Ta exp{-a;7(2a;i)}, 
,{x) = ni{l + {x/xi)'^) -I- n2 (l + (x/x2f) 



(9) 
(10) 

where xi = 1.54 x 10^^ cm = 5 kpc, X2 = 7.19 x 10^^ cm = 
23.3 kpc, ai = 1.518, a2 = 0.705, Ti = 2.78 x lO'^ K, and 
Tz = 8.997 x 10^ K, m = 0.089 cm-^ na = 0.019 cm-^ 
(see Fig. [5}, we compute the gravitational acceleration due 
to the cluster's gravitational potential from the assumption 
of hydrostatic equilibrium (HSE) , 

IdP 

--r- = 9{x), 
p ax 

where P — nksT is the thermal pressure, and g{x) is the 
gravitational acceleration. 

At the beginning of the simulations we set the temper- 
ature field to be uniform, Tq = 3 x 10^ K, and the initial 
density distribution is determined from the requirement of 
HSE. Since HSE fixes only the gradient of the density, we are 
left to choose the peak density, po, ad hoc. From a number 
of one dimensional test simulations (without AGN heating) 



we determined the optimal value as po = 5.15 x 10 g cm~ , 
which provides the best fit to the observational profiles on 
the outskirts of the cluster at time t ~ 10^" yr a fter the be- 
ginni ng of the simulation (see also simulations bv lPope et al.l 
l2005h . 



2.3 Introduction of Bubbles 

We use a single fluid in flash to represent both ambient 
ICM and the gas inside bubbles, unlike the two-fl uid ap- 
proach used by Briiggcn ct al. (200^ andlRuszkowsk Fet al.l 
(200i). The bubbles are introduced into the simulation when 
the temperature of a computational cell in the vicinity of 
the centre falls below T — 1.5 keV. This trigger temper- 
ature serves as a simple feedback mechanism that relates 
cooling flow parameters to the heating response. The value 
1.5 keV was selected to be close to the observed temperature 
minimum, see Fig. [2] 

The temperature field perturbations associated with the 
bubbles are introduced using a modified version of the per- 
turbation module util/perturb, which is supplied with 
FLASH. We position the centres of the perturbations sym- 
metrically on opposite sides of the centre of the grid (clus- 
ter), at distances of xq = 1.5ro — 1.5 x lO'^^ cm = 4.86 kpc, 
with randomly specified orientation angles of the line con- 
necting the bubbles and the cluster centre. The locations for 
the bubbles are refined up to the maximum refinement level 
prior to the introduction of the bubbles. 

The temperature of all computational cells with coor- 
dinates b, satisfying the inequality ro < |b — Ci|, where Ci 
{i — 1,2) are the position- vectors of the bubble centres, 
is modified using the following algorithm. Each computa- 
tional cell b is sub-divided into 4x4x4 sub-cells (bs, 
s — 1, . . . ,64). The temperature of each sub-cell is set to 
Tb = 5 X 10^° K if ro < |bs — Ci|, and is left unchanged 
otherwise. The temperature of a computational cell is then 
calculated as a volume average of temperatures of its sub- 
cells, r(b) = l/64 5^f T(b,). This algorithm provided a 
sharp transition boundary between the bubble and the ICM 
(so-called top-hat perturbation profile), but is sufficiently 
robust not to cause any numerical difficulties. 

To make the temperature perturbations thermodynam- 
ically consistent, we choose a pressure profile, and calculate 
the corresponding densities using the equation of state for 
an ideal gas. We tested two profiles: the pressure inside the 
bubble matches the pressure of the ambient medium {i.e. 
the pressure field is not changed by the perturbation), and 
a constant pressure profile, where the value of the constant 
pressure is calculated as a volume average of the pressure 
field inside the perturbation region before the introduction 
of the bubbles. However, we find that since the sound speed 
inside the bubble is much higher than the sound speed of the 
ambient medium, in the former case the pressure inside the 
bubble quickly reaches a constant profile, sending a weak 
sound wave into the ICM. In order to avoid this additional 
artificial disturbance we used constant pressure profiles in 
our setup. Further analysis (see § |3]) showed that the dif- 
ference in the initial conditions does not lead to a serious 
discrepancy in the properties of the bubbles. The minor ef- 
fects caused by it will be discussed in § [l] As a first order 
approximation, however, the difference between the param- 
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Figure 2. Left panel: observed temperature and density profiles of the Virgo cluster from G04. The lines are the analytical fit, see 
equations ^ and JTOj. Right panel: gravitational acceleration as a function of distance from the centre, derived from HSE conditions 
using the best-fit lines. 



eters of the bubbles in these two cases can be viewed as an 
intrinsic scatter caused by our numerical scheme. 



3 SIMULATION RESULTS 

In this section we present an analysis of the data from sim- 
ulations with two kinds of the initial pressure profiles for 
bubbles, model A has the initial pressure profile matching 
the ambient pressure and C has the initial constant pressure 
profile. 

The core of the cluster met the criterion for injection 
of bubbles after t = 706 Myr, so the bubbles are introduced 
into an established cooling fiow. At this stage the ambient 
temperature of the ICM at the distance xq = 1.5ro from 
the centre is T{xo) « 2.62 x lO'' K, density p{xo) « 2.83 x 
lO"^*^ g cm"^, and the cooling time At « 1.7 Gyr. Due 
to the large cooling time we can ignore the effects of the 
radiative cooling on the parameters of the bubbles. 

In our purely hydrodynamical simulations we need to 
resolve the smallest scales of the instabilities in order to ac- 
count properly for the small scale mixing of the hot gas with 
the ambient ICM. The key parameter in the stabilisation of 
the bubble/ambient medium interface against the RT in- 
stability is the surface tension, which can be provided by 
the magnetic field lines parallel to the sur face of the bub- 
ble. As it was shown b y lKaiser et al.l l|2005h (hereafter K05) 
and iDe Youngj l|2003h even weak magnetic surface tension 
(magnetic fields of ~ O.l/iG) stabilise the bubble surface on 
scales up to 0.1 kpc at the distance of 5 kpc from the cen- 
tre of the cluster (see Fig. 3, K05). By limiting the smallest 
cell size to h = 0.16 kpc we essentially assume that the gas 
will be fully mixed on the scale ~ /i on the time scale of 
h/cs ~ 10^'' s, which is approximately equal to the growth 



time of the RT instability of the scale h in the absence of vis- 
cosity (see Fig. 2, K05). Although a fully consistent model 
should include magnetic fields to account for the stabilisa- 
tion of the bubble surface (Ruszkowski ct al. 2007), these 
simple estimates show that the numerical resolution of our 
model is consistent with the expected dynamics of the sys- 
tem, as follows from the analytical stability analysis, and 
numerical simulations with 'random' magnetic fields. 

3.1 Bubble Surface 

As the bubbles ascent in the ICM their shape and position 
changes. In order to study their parameters we have used 
a simple technique for identification of the material inside 
bubbles at any stage during their evolution. Firstly, we note 
that the averaged radial entropy index (cr = p~^^^T) pro- 
file for the ICM, can be well fitted with a linear function, 
(a) — a^x + (Jb, where a; is a distance from the centre, at 
any time during the simulations in both models. Secondly, 
by constructing the entropy profile using the means over 
large enough spherical shells (i.e., using low resolution pro- 
filing) we can abate the contribution from the high entropy 
regions associated with the bubbles (Vbubbic/Vishoii ^ !)• 
The variances of the values of a in each shell are used as 
an error estimate for the values of the means in the fitting 
routine, which further reduces infiuence of the high entropy 
bubbles on the result of the linear fitting. We mark as "bub- 
ble" any computational cell with an entropy index above 
1.5(aaa:: -|- (Tb), where x is the distance of the cell to the cen- 
tre of the cluster (grid), and aa,b are the fitted coefficients 
for the current time. 

The surface of the bubbles becomes irregular due to RT 
instabilities, and small bits of the bubble are shredded away 
as it rises through the ICM and fragments. In such circum- 
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stances a measure of the size and position of the bubble 
material has to be essentially statistical, and not just purely 
geometrical. 

To quantify the size of the bubble we measure its ex- 
tension in the direction of ascent, _R||, and its size in the 
direction perpendicular to the direction of ascent, R±. For 
a purely spherical bubble 7?|| = R±. When the bubble has 
the shape of a torus we measure both external and internal 
radii of the torus, -R±(max,min) = max, min(_Rx), using the 
following algorithm. The position vector of a "bubble" cell, 
Ri, is split into the component aligned with the direction 
specified by a vector connecting the centre of the cluster 
and the initial position of the bubble, R||i (the direction of 
ascent), and the vector perpendicular to the first one, going 
through the centre of the cell, Rxi, so that, Ri = R||i+Rxi. 
We then compute histograms (using 100 bins) of |Rxi| and 
|R||i|, using the values of the corresponding vectors R(||_x)i 
for all cells of a single bubble. We use the cut off at the 
level of 0.1 of the maxima of the histogram to find the ap- 
proximate extent of the bubble material (and exclude any 
trailing and overtaking bits; we select the bins below the 
cut-off which are closest to the peak of the histogram), that 
gives the inner and outer radii of the torus, -Rxmax,min, or, 
in the case of the 'quasi-spherical' bubble, two measures of 
the radius, _Rx = -Rxmax, and 7?|| = 0.5(i?||max - -R||min)- 
The means of the |-R(||,x)i| with the exclusion of the values 
below the cut-off, correspond to the distance of the bubble 
from the centre of the cluster (in the case of |R||i|), and a 
mean radius of the torus (in the case of |Rxi|). 



3.2 Morphology 

The surfaces of bubbles rising in the ICM corrugates due 
to Reilegh- Taylor (RT) instability, however, this does not 
result in an immediate dispersion of the bubbles, see Fig. [T] 
In fact, the mass diffusion due to the RT instability and 
numerical diffusion remains quite small and the bubbles in 
simulations A and C remain traceable for a long time. We can 
still find traces of the old bubble material after the second 
pair of bubbles was injected (> 100 Myr). 

We find that during the ascent the size of the bubbles 
does not change adiabatically according with the change of 
the ambient pressure with distance from the centre of the 
cluster, 

Vix) = Vo{P{x)/Po)-'/'', (11) 

where P{x) is the pressure profile, with 7 — 5/3. Also the 
velocity of the bubbles does not reach a terminal level, as it 
is us ually expected from a simple theoretical analysis (see, 
e.o.. lEnfilin fc Heinz|[200^ . 

As shown in Fig. [S] the volume of the bubbles keeps 
growing after the injection and then it plateaus. The growth 
of the volume of the bubble (which is identified as the re- 
gion of high entropy) is the consequence of several different 
processes. One of them is the adiabatic expansion due to the 
pressure change with distance from the centre of the cluster 
as given by equation (|lip . another one is the entrainment 
of ambient plasma through the surface of the bubble (mix- 
ing), and the third one is the enlargement of the bubble 
due to the hydrodynamical Kutta-Zhukovsky forces as ex- 
plained below. Full quantitative description of the evolution 
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Figure 3. Change of the volume of bubbles with height. Note 
that due to discretisation (finite resolution of the grid) and the 
algorithm we use for identification of the bubble cells during the 
analysis, the initial volume of the bubble is less than the geometri- 
cal value Vb = 4/37rrQ. The dash-dotted line shows the predicted 
isotropic adiabatic expansion of the spherical bubble according to 
equation JTTJ. The initial volume is taken to be equal to the initial 
volume of bubble determined by the bubble finding algorithm. 

of the bubble parameters is out of scope of this work, an d is 
addressed in the follow-up article l|Pavlovski et al.ll2007l ). 

In our numerical experiment the entrainment of the 
ICM through the boundary of the bubble is a consequence 
of the instabilities that develop on the surface of the bubbles 
(see Fig. [7J , and, unavoidably, the numerical diffusion. In re- 
ality the diffusion of the plasma across the boundary of the 
bubble is controlled by the magnetic field, and the resulting 
change of the volume of the bubble is likely to be differ- 
ent both from the prediction given by equation (|lip . and 
the results of our numerical model. It is important to note, 
however, that any amount of mixing will lead to violation 
of the adiabatic assumption. Unless the mixing process is 
understood in greater detail, and other adiabatic processes 
like Kutta-Zhukovsky forces are taken into account, the es- 
timate of the AGN energy input from the estimates of the 
total volume of the bubbles is likely to be inaccurate. 

The initial growth of the volume of the bubbles corre- 
sponds to the deceleration phase in their ascent, which fol- 
lows the very sharp initial acceleration, see Fig. The slight 
re-acceleration of the bubbles at x ~ 3ro corresponds to the 
change of the morphology of the bubbles from spherical to 
toroidal, after which the bubbles continue to decelerate. As 
explained below, the velocity of the bubble is tightly linked 
with the increase in its volume by a simple physical mecha- 
nism. 

The initial expansion of the bubble is not isotropic, 
see Figs. [5] and [B] The volume of a bubble increases due 
to the sideways expansion (expansion in the direction per- 
pendicular to the direction of ascent, _L-direction) , whereas 
the size of the bubbles in the direction of the ascent (||- 
direction) stays roughly constant. The ellipticity of the bub- 
bles, e = ^1 — (ii||/iix)^, quickly reaches saturation at a 
level of e « 0.9 in both cases. 

The flattening of the bubbles can not be attributed to 
the viscosity, since no physical viscosity is modelled. In pres- 
ence of viscous stresses they would act to squeeze the bubble 
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Figure 4. Change of the ascent velocity of the bubbles with 
height. 
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Figure 5. Change of the size of the bubbles in the direction of 
the ascent (| |-direction), and the ellipticity, ^1 — \/ Ri_ of 
the bubbles. 

in the | |-direction. They would resuh in an increase of pres- 
sure inside the bubble, which would act to restore the equi- 
librium. This process can result in oscillations of the surface 
of the bubble (by analogy with oscillations of the buoyant 
bubbles in a lava lamp), but is unlikely to lead to any net 
increase in their volume. 



3.3 Kutta-Zhukovsky forces 

The structure of the flow through the bubble has the ge- 
ometry of a vortex ring right from the very early stages, 
when the bubble is still roughly spherical, see Fig. [S] This 
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Figure 6. Change of the size of the bubbles in the direction 
perpendicular to the direction of ascent (±-direction). 




Figure 7. Surfaces of the bubbles. The bubbles closest to the 
centre correspond to the time of the injection and are spheri- 
cally symmetric. The later stages correspond to times 18 Myr 
and 98 Myr after the injection for (model c). Note the 'finger- 
like' corrugation of the surfaces of the bubbles. 

circumstance has a profound effect on the evolution of the 
bubble. 

As a hot bubble starts to ascend, it can retain its spher- 
ical shape only if the fluid motion inside the bubble corre- 
spond to a perfect fluid dipole^ (by analogy with so-called 

^ There is a large amount of literature about (quasi-) two- 
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Figure 8. The flow around the spherical bubble shows the struc- 
ture of a vortex ring. The visualisation shows part of the surface 
of the bubble (the upper half is cut off) in simulation C, 2.3 Myr 
after the injection. The arrows show the direction of the velocity 
in a thin slab just below the cut plane. The tubes trace the stream 
lines close to the centre of the right vortex on the cut plane. 



Hill's solution, see, e.q., ISaffman|[l995h . If during the as- 
cent the spherical symmetry breaks down (i.e., R± becomes 
larger than R\\), the fluid motion inside the bubble corre- 
sponds to the circulation around a stretched dipole - the 
vortex ring. 

In the beginning we can consider the inner radius of 
the vortex ring (torus) to be equal to zero (the bubble is 
roughly spherical), and the velocity of ascent to be equal 
to the maximum ascent velocity, see Fig. O Due to drag the 
speed of the ascent and of the vortex flow will decrease. Since 
the vortex moves with respect to the ambient medium each 
element of the vortex is subject to the Kutta-Zhukovsky^ 
force (for an introductio n to aerodynamic forces see, e.g., 
iLandau fc Lifshit3 (|l987l ) § 38, and appendix|3, which acts 
to expand the vortex radially. This force acts to expand 
the bubble in the _L-direction, eventually creating a hole 
in the middle of the bubble - the bubble becomes torus 
shaped. As the vortex ring expands radially, an additional 
component of the Kutta-Zhukovsky force pushes the bubble 
down, reducing its ascent velocity. Indeed, with rapid vortex 
expansion the total velocity vector of the elements of the 
vortex ring will no longer be aligned in the ||-direction, but 
will have a component in the _L-direction as well, see Fig. [Q] 
This results in the Kutta-Zhukovsky force pointing down 
in the ||-direction, in the direction opposite to the velocity 
of ascent. This force will act to reduce the velocity of the 
ascent, i.e., it counteracts buoyancy. 

Let R — {R\\) be the radius of the torus, and r the 
radius of its cross-section, see Fig. |§] If a is the central an- 
gle, so that Tida is the length of an element of the vortex 
ring, then this element will be subject to a Kutta-Zhukovsky 



dimen s ional vortex dynamics, s ee, e .g., iTurner] lll957l): [Mortonl 
lll960l '): lTurneil lll969l '): lFraenke]| l ll972l) . see also lAfanasvevI \200^ 
for a review of experiments. 

^ Zhukovsky's surname is sometimes spelt as Joukovsky or 
Joukowsky in the literature. See, e.g., Joukowsky transform, also 
Kutta-Schukowski, Kutta-Joukowski, etc. 




Figure 9. The forces acting on an element of the vortex ring: 
_Fb - buoyancy, - drag, F]^^ - Kutta-Zhukovsky force. Dashed 
arrows show the components of the velocity: the ascent velocity 
U = X and the expansion velocity u = R. 



force, 



1/2 



Rda, 



(12) 



where 7? — dR/dt is the speed of the vortex expansion, F 
is the vortex circulation (see appendix , pa is the density 
of the ICM, and U' — U ~ v' , U is the ascent vel ocity, v' is 
the s elf-induced (vertical) velocity of the vortex (|Batchelorl 



F , R 

V ~ -; log ■ 



The Kutta-Zhukovsky force points at an angle (3 to the _L- 
direction. 



sin (3 ■ 



cos P 



R/U' 



l + [RIU'y 
1 



(14) 



1 + 



So, the vertical component of the Kutta-Zhukovsky force, 
acting on the whole vortex ring is, 

= 27rEr7?pa. (15) 

In Fig. [To] we plot the ratio of the vertical component 
of the Kutta-Zhukovsky force to the buoyant force acting on 
the bubbles as a function of their distance from the cluster 
centre. To estimate the circulation of the vortex we calcu- 
lated the mass-weighted mean of the velocity of the bub- 
ble material in the rest frame of the bubble. This is clearly 
just an estimate, and the resulting plot has some scatter 
as a result. It does, however, illustrate the main points dis- 
cussed above, and shows that the magnitude of the Kutta- 
Zhukovsky force is comparable to the buoyant force. A quan- 
titative model of the ev olution of the bubble s is out of the 
scope of this work (see iPavlovski et al. I l2007l, for more de- 
tails) . 



8 Georgi Pavlovski et al. 



1 .0 




dis;ance / r„ 



Figure 10. Change of the ratio of the magnitude of the ||- 
component of the Kutta-Zhukovsky force to the magnitude of 
the buoyant force. 

4 DISCUSSION 

In the present work we have omitted the phase during which 
the jets inflate the bubbles, which hmits apphcabihty of this 
model to the interpretation of observations. However, the 
dynamics of the flow around bubbles created by a jet is 
likely to be similar, once the bubbles become buoyant. Tak- 
ing the jet injection into account it is quite possible that 
the vortex ring structure will form much earlier (during the 
injection) and the flow around the bubble will have a larger 
initial value of circulation, F, which can be induced by the 
viscous momentum transfer of the kinetic energy of the jet, 
as suggested by the laborator y experiments of small scale 
jets (see, e.g.. I Afanasvev|[2006l ) . Kutta-Zhukovsky forces are 
then likely to be even more important for the dynamics and 
the overall energy balance. 

The parameters of the bubbles with different initial 
pressure proflles (models A and c) do not differ greatly, see 
Figs. [3] - in] The only noticeable difference is that the bub- 
bles in simulation A develop a comparatively large secondary 
torus. This secondary vortex ring (torus) is shredded from 
the main vortex soon after the bubbles morph into tori, it 
separates from the much bigger main torus, and travels in 
front of it. The volume of this secondary torus is small com- 
pared to the volume of the bubble, and its contribution to 
the size of the bubble was discarded (it has generally fallen 
below the 10 per cent cut-off, see § I3.1|l during the calcu- 
lations. The secondary vortex ring quickly mixes with the 
ambient ICM and dissipates. The fact that this vortex ring 
has a larger volume in case A is reflected in the evolution 
of the total volume of the bubbles, see Fig. [3] The volume 
of the A-bubbles is getting smaller than the volume of the 
C-bubbles once their secondary tori mix with the ICM. It 
also results in a slightly smaller i?x size of the A-bubbles at 
later stages of the evolution, see Fig. |6] 



4.1 The role of missing physics 

High temperature plasma in the absence of the magnetic 
fields is thou ght to be thermally conductive and viscous 
l|Spitzeij|l963 ). The diffusivity coefficients are sensitive func- 



tions of temperature, cx T^'^ . This ensures that any tem- 
perature gradient on scales larger than the mean free path 
of an electron is smoothed out very efficiently by the con- 
duction, and fluid instabilities are damped by the viscosity 
iKaiser et al.ll2005l ). 

We have attempted to simulate the dynamics of bub- 
bles in a viscous and thermally conductive ICM. Our exper- 
iments showed, however, that unless thermal conduction is 
very efficiently suppressed on the boundary of the bubble, it 
leads to a very fast erosion of the bubble material (although 
no fluid instabilities develop). This is easy to understand 
from the following simple estimate. Given that the amount 
of heat in the bubble is of the order AQ ~ pV « 10^^ erg, 
the difference in the temperatures of the plasma inside the 
bubble and the ambient ICM AT ~ T 10^ K, the coefli- 
cient of the thermal conductivity for plasma at T = 10^ K 
is K ~ 10^^ erg s~^ cm~^ K~^, then the timescale of the 
evaporation of the bubble is At = AQAI/{k,ATS) ~ 1 Myr, 
where we let AZ ~ i?, and S ^ R'^. Therefore it is not at all 
surprising that we found that the bubbles were fully mixed 
with ICM at a distance of 3 . . . Aro from the centre of the 
cluster. Without the inclusion of magnetic fields into the 
numerical model it is impossible to draw a definitive conclu- 
sion about the evolutio n of the bubbles in a diff usive ICM. 
Magnetic field draping (|Ruszkowski et al]l2007l ) is likely to 
play a significant role in the dispersion of the bubbles in 
thermally conductive ICM. 

We note, however, that the above discussion of the KZ 
forces does still remain valid even in our 'superdiffusive' 
cases (see appendix iBl for further details), i.e. the bubbles 
still morph into tori, and the vortex motion is present. 

The morphology of the observed bubbles in the Perseus 
cluster w as used as an indicato r of a possible highly viscous 
ICM in (|Revnolds et al.ll2005h . We would argue, however, 
that the flatten i ng of the bubble, which was explained by 
[Reynold's et all l|2005h as a result of viscosity, has an ad- 
ditional explanation. It is the Kutta-Zhukovsky forces, not 
viscosity, that determine the radial size of the bubbles, and 
it is impossible to infer the importance of viscous effects in 
the ICM from the morphology of the bubbles alone. Viscos- 
ity acts to suppress RT instabilities but does not prevent 
morphing of the initially spherical bubbles into tori. At this 
stage the presence of a significant viscosity in galaxy clusters 
should be regarded as speculative, unless directly confirmed 
by observations. 



5 CONCLUSIONS 

In the present work we have used a numerical scheme based 
on the FLASH code to produce three dimensional simulations 
of AGN bubble heating of the ICM, with an emphasis on the 
dynamics of the bubbles. Our setup includes a sophisticated 
cooling function to help establish the initial cooling fiow, fol- 
lowed by the introduction of the high contrast AGN-bubbles, 
which dynamics we have analysed. 

The data presented in this article are consistent with 
the results of previous num erical simulations of the AGN 
bubb l e heating of the ICM llGardiiiil l2007l : [Reynolds etall 
l2005l : iBriiggen fc Kaisej|200^ However, our analysis has 
indicated additional important phenomena linked to the dy- 
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namics of the gas of the bubbles in the ICM. Our main re- 
sults are as follows. 

(i) Bubbles change shape and transform into tori not be- 
cause of hydrodynamical instabilities, but due to the radial 
Kutta-Zhukovsky force. 

(ii) Bubbles do not reach a terminal velocity defined by 
the balance between the drag force and the buoyancy. The 
velocity of their ascent is also affected by the the vertical 
Kutta-Zhukovsky force and the self-induced velocity of the 
vortex. 

(iii) Viscosity does play a role in stabilising the surface of 
the bubbles against instabilities, but it does not alter their 
overall 3D (torus-like) morphology. 

(iv) The volume of bubbles can lead to an erroneous esti- 
mation of the energy supplied by the AGN when using equa- 
tion (jlip . The expansion of the bubbles is not only due the 
change of the ambient pressure but also due to the Kutta- 
Zhukovsky forces, and non-adiabatic turbulent mixing. A 
more accurate model of plasma (including magnetic fields) 
is needed in order to make realistic prediction about the 
entrainment of the ICM into the bubbles. 

(v) Thermal conduction must be highly suppressed on the 
boundary of the bubbles, otherwise it leads to a rapid erosion 
(mixing) of the bubble material. 

A better knowledge of the physical state of the ICM is of 
paramount importance for further progress in understanding 
the interaction of AGN outflows with the ICM. 
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APPENDIX A: NATURE OF THE 
KUTTA-ZHUKOVSKY FORCES 

Circulation is the line integral around a closed curve of the 
fiuid velocity. If v is the fiuid velocity and ds is a unit vector 
along the closed curve C, the circulation is given by, 

T = ^ V • ds. 

The units of circulation are [cm^ s~^]. The Kutta-Zhukovsky 
theorem states that the lift force acting per unit span on a 
body in an inviscid fiow field can be expressed as the product 
of the circulation about the body, F, the fluid density, p, and 
the speed of the body relative to the free-stream, U, 

f = UpT. 

This equation applies both around airfoils, where the cir- 
culation is generated by airfoil action, and around spinning 
objects, experiencing the Magnus effect, where the circula- 
tion is induced mechanically. The latter case is governed by 



the same basic principles as apply to the travelling vortices. 
A spinning object {e.g., a cylinder rotating around its axis 
of symmetry, or an eddy) creates a motion of the fiuid in the 
boundary layer around it {e.g., via viscous drag) - the cir- 
culation. If such a rotating object travels through the fluid 
{e.g., the rotating cylinder travels through the fiuid at rest 
in the laboratory frame of reference, with its axis of sym- 
metry remaining parallel to itself), on one side of it, the 
velocity of the fiuid in the boundary layer will be in the 
same direction as the velocity of the surrounding fluid near 
it. On this side the resulting velocity of the fiuid will be 
larger than the travel velocity. On the other side of the ob- 
ject, the velocity of the fluid in the boundary layer will be 
in the opposite direction of the velocity of the fiuid near it, 
and the resulting velocity of the fluid will be smaller than 
the travel velocity. According to the Bernoulli's theorem the 
quantity p + v'^ /2 = const along fiow lines, so the pressure 
p will be lower on one side of the rotating cylinder than on 
the other causing an unbalanced force at a right angle to 
the travel velocity. 



APPENDIX B: NUMERICAL EXPERIMENTS 
WITH THERMAL CONDUCTIVITY 

In order to account for the suppression of the diffusive pro- 
cesses in the presence of chaotic magnetic fields the diffusive 
coefficients for conductivity, Ks, and viscosity, r/s, are reduced 
by a certain factor, < /s < 1. In the case of AGN blown 
bubbles this suppression is likely to be much stronger inside 
the bubbles (due to the associated strong magnetic field), 
resulting in the suppression coefficient being not spatially 
uniform, 

K = /s(x)ks, (B1) 
V = fs{^)Vs, (B2) 

/.(x) = |-^=^' «-tside bubbles, ^^^^ 
l/s2, inside bubbles. 

The factor fsi is a simulation parameter^. We have tested 
a range of possible values: fsi = [0.1,0.3,0.6]. In order to 
identify the computational cells that are located inside the 
bubbles and to apply a different suppression factor to them, 
we probe cells inside the spherical central region of the size 
X7 = 7ro = 7 X lO'^^ cm = 22.7 kpc, where ro is the initial 
radius of the bubbles, and X7 — 7ro is an estimate of the 
maximum height the bubbles are likely to reach, starting 
at the distance of 1.5ro from the centre of the cluster (see 
also ij |2.3p . At each time-step we compare the entropy index, 
a = Tp~^^^ , of the cells inside the |x| < X7 region with the 
cutoff value, ar. The cutoff value is computed for each time- 
step as, CT7 — 1.5((T)2-.^<|x|<a;j+Aa;, which Is 1.5 times the 

^ Note that the suppression factors for viscosity and ther- 
mal conductivity are likely to be similar, despite the fact that 
the corresponding Larmor radii are different by several orders 
of magnitude. In the tangled magnetic field with the coher- 
ence scale the Rechester- Rozenbluth distance (for details see 
iNaravan fc MedvedevI I2OOII) for electrons is given by, Z/j^j^ ~ 
ln(ZB/?'e) ~ IOZb, where is Larmor radius for electrons. For 
ions this distance is very similar Lrr ^ 25iBi implying similar 
suppression coefficient for both transport processes. 
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Figure Bl. Schematic of the growth of a conductive layer around 
the bubble due to the discretization of the temperature gradient. 
The vertical lines represent boundaries of the computational cells. 
The solid line represents the temperature profile at the edge of the 
bubble, where Tj, is the temperature of the gas inside the bubble, 
and Ta is the temperature of the ambient medium. The dashed 
line represents the corresponding value of the thermal conduction. 
The suppression coefficient /si is applied to the gas at tempera- 
ture Ta, and /s2 to the gas at temperature Tj, . The panel on the 
left represents the initial profiles. Thermal conduction works to 
reduce the temperature gradient, and creates a number of cells 
with an intermediate values of temperatures (entropies) and high 
conductivities (right panel). 




Figure B2. Growth of a conductive layer around the bubbles in 
a 3D simulation. The visualisation shows surface cuts through the 
scalar field equal to the logarithm of the thermal diffusivity. The 
resulting surfaces are colour-coded and warped according to the 
value of the scalar. Semi-transparent surfaces show the edge of the 
bubbles. The data are taken from the simulation with /gi = 0.1: 
on the moment of injection (left), 8 Myr later (right). The size of 
the box is given in kiloparsecs. 



average of the entropy index in a spherical shell with radius 



x-j, and thickness As 
then define fs2 as, 



0.5ro = 5 X 10^' cm = 1.6 kpc. We 



/s2 = 



X 10-^" (2p 

* Arr 



if (J > a?, 
otherwise, 



(B4) 



where the factor (2(77/3(t)^ '' ~ (Tk/Tb)^'^ accounts for the 
increase of the conductivity and viscosity coefficients due to 
the large temperature of the bubbles, Tb, compared to the 
temperature of the ambient gas, Ta. This scheme is designed 
to effectively suppress diffusive process inside the bubbles. It 
could not, however, prevent diffusion of the bubble material 
and its mixing with the ICM. 

In our diffusive simulations discretization of the temper- 
ature gradient leads to formation of a layer of thermally con- 
ductive plasma around the bubble, as illustrated in Figs. IBll 
and lB2l This layer is essentially an artifact of our numerical 
scheme for the suppression of the diffusive processes, and is 
not reflective of the real physical properties of the ICM. 

The problem with discretization is going to be present 
in the two-fluid approach as well (by two-fluid approach we 



mean separate fl uids for the bubble and the ambient medium 
as used by, e.g.^ iBriiggen et ai] l|2005l )). We believe it will 
lead to a qualitatively similar behaviour, since the suppres- 
sion coefficient for the diffusion in this case is also based on 
an ad hoc cut off based on the concentration of the bubble 
fluid in a computational cell. 

Without a more advanced model, which needs to in- 
clude magnetic fields, and self-consistently accounts for the 
suppression of the diffusive processes, it is impossible to say 
whether or not this is a realistic representation of the real 
state of the ICM. 

We find that the dynamics and morphology of the 
bubbles in our diffusive simulations also governed by KZ 
forces, since the circulative flow around the bubble is always 
present. 
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